Liquid-vapor transition from a microscopic theory: Beyond the Maxwell construction 
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A smooth cut-ofT formulation of the Hierarchical Reference Theory (HRT) is developed and applied 
to a Yukawa fluid. The HRT equations are derived and numerically solved leading to: the expected 
renormalization group structure in the critical region, non classical critical exponents and scaling 
laws, a convex free energy in the whole phase diagram (including the two-phase region), finite 
compressibility at coexistence, together with a fully satisfactory comparison with available numerical 
simulations. This theory, which also guarantees the correct short range behavior of two body 
correlations, represents a major improvement over the existing liquid state theories. 
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Given that the physics of the hquid- vapor phase tran- 
sition is a textbook topic in thermodynamics, it is dis- 
appointing that so far no liquid-state theory has been 
able to describe it satisfactorily. In fact, when faced with 
phase coexistence, mean-field approaches such as (gen- 
eralized) van der Waals theories give a non-convex free 
energy, while integral equations [l| fail to converge al- 
together in some domain inside the coexistence region. 
In both cases, the coexistence boundary must be recov- 
ered via the Maxwell construction, namely, by imposing 
thermodynamic equilibrium a posteriori. Aside of first- 
principles considerations, it should also be noted that 
this procedure may turn out to be cumbersome to im- 
plement, especially in the case of integral equations, be- 
cause of the need to circumvent the forbidden domain 
and the ambiguities entailed by this procedure for the- 
ories lacking thermodynamic consistency. When dealing 
with mixtures of fluids, this becomes a serious hindrance 
to the theoretical determination of the phase diagram. 

The approach which comes closest to a realistic de- 
scription of the liquid- vapor transition is the Hierarchical 
Reference Theory of fluids (HRT) [1], a genuine micro- 
scopic approach which implements the Renormalization 
Group (RG) method within the liquid state framework. 
HRT is able to reproduce non classical critical exponents, 
scaling laws and, what is most relevant for this discussion, 
a rigorously convex free energy, so that flat isotherms at 
coexistence naturally emerge from the theory, leading to 
accurate predictions of phase boundaries. However, the 
implementation of HRT developed so far is plagued by an 
unphysical divergence of the compressibility on the coex- 
istence boundary In liquid-state theory jargon, HRT 
forces the spinodal and the binodal curve to coincide. If 
flat isotherms are obtained at this price, one may have 
some doubts about how reliable the description of the 
first-order phase transition given by HRT really is. This 
also prevents one from achieving any information on the 
possibile occurrence of metastable pure phases which, in 
the usual approaches, are located between the binodal 
and the spinodal. 

In this Letter we show that a novel implementation of 



HRT based on a smooth cut-off procedure [J] , is able to 
reproduce correctly the physics at the phase boundary 
and provide a criterion to discriminate between unstable 
and metastable states, while also improving on the repre- 
sentation of short range correlations. This qualifies HRT 
as the only liquid state theory able to provide a satisfac- 
tory description of the liquid-vapor phase transition. 

Following the lesson of the RG [11, HRT is based on 
the gradual introduction of (density) fluctuations start- 
ing from short wavelengths. The two body potential 
v{r) is first split into the sum of a (repulsive) refer- 
ence part viiir) and a (mostly attractive) tail w(r). The 
properties of the reference system, usually a hard sphere 
fluid, are assumed to be known, and a sequence of inter- 
mediate systems is introduced, labeled by a parameter 
t G (0, oo) "interpolating" between the reference and the 
physical model. The interaction potential of the i-system 
is vt{r) = vji{r) -\- Wt(r) where the Fourier components 
Wt{k) of Wt(r) are strongly suppressed for wave vectors 
k < e~*. The change in the free energy of the system 
when the parameter t undergoes an infinitesimal change 
can be evaluated exactly, leading to the HRT differential 
equation. 

The purpose of the procedure outlined above is to sup- 
press the liquid- vapor transition throughout the whole se- 
quence of intermediate t-systems, due to the long range 
repulsive tail present in Wt{r). Like in the RG, the long- 
wavelength fluctuations which drive the phase transition 
must be allowed to develop only in the t —>■ oo limit, 
when Wt{r) tends to the physical attractive interaction 
w{r). The sharp cut-off formulation of HRT, which has 
been successfully applied to several model systems [6|, 
fulfills this condition. The sharp cut-off is defined by the 
choice Wt{k) = w{k) for k > e~* while Wt(fc) — for 
k < e~* and gives rise to an HRT equation which, close 
to the critical point and at long wavelengths, reproduces 
the Wegner-Houghton RG scheme 0] . As previously dis- 
cussed, this formulation suffers from some deficiency close 
to the first-order boundary, where a diverging compress- 
ibility is predicted by the sharp cut-off HRT. Remark- 
ably, it has been recently shown, in the framework of 
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scalar field theories, that this unphysical behavior can 
be eliminated by changing the cut-ofF procedure, i.e. by 
modifying the definition of the intermediate potentials 
Wt (f) [1, Q . If this property survives in the framework 
of a liquid state theory is not known and it is addressed 
here. In this Letter, we adopt the following form: 



Wt{r) ~ w{r) 



■ilj{t)w{re *) 



(1) 



where c? = 3 is the space dimensionality and ipit) is a 
decreasing function of t with ?A(0) = 1 and asymptotic 
behavior ip{t) oc e~^* for t — > oo. The precise definition 
of ^p{t) will be discussed later. Note the t dependence of 
the range of the second (weakly repulsive) contribution 
to wt{r): in the large t limit, the amplitude of this term 
decreases while its range grows. 

The differential equation expressing the change in the 
free energy when the parameter t (and hence the interac- 
tion vt) is changed follows from first-order perturbation 
theory 0: 



dt 2 ^ ^ ' dt 



(2) 



where At is (— /?) times the Helmholtz free energy per 
unit volume of the fluid interacting via Vt{r) and gt{r) 
is the corresponding radial distribution function, while 
(j)t{r) = —Pwt{r). This equation is formally exact but, as 
usual in the HRT approach, it requires a closure relation 
expressing the two body correlations in terms of the free 
energy. Analogously to the standard implementation of 
the sharp cut-off HRT Q , we adopt the following Mean 
Spherical Approximation (MSA)-like representation [l^] : 

gt{r) = for r < 1 (3) 

ct{r) = (j)t{r) + \t (l>{r) for r>l (4) 

where the direct correlation function Ct{r) is related to 
gt{r) by the usual Ornstein-Zernike (OZ) equation [lo| . 
The key feature of this closure is the presence of the A* 
parameter which is introduced in order to enforce the 
compressibility sum rule: 
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dr ct{r) 



(5) 



Equations ^SSij form, together with the definition ([T]), a 
closed set of integro-differential equations for the ther- 
modynamics and correlations of the model. The numer- 
ical solution of this problem poses a highly demanding 
computational tasks, which may be however simplified 
by choosing a particularly favorable model system: the 
Yukawa fiuid, i.e a hard sphere fluid with attractive tail 
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where T is the dimensionless temperature and z is the in- 
verse range. Lengths are normalized to the hard sphere 



diameter. This one-parameter family of interaction po- 
tentials actually represents one of the most studied sys- 
tems in liquid state theory: due to the simple analytical 
form oiw{r) and the flexibility due to the tunable inverse 
range parameter z, it provides a reasonable description 
of simple fluids (for the celebrated choice z = 1.8) as well 
as colloidal suspensions, where typically z ^ 1. A ma- 
jor advantage of this particular form of the interaction, 
follows from the availability of an exact analytical solu- 
tion of the OZ integral equation with the ansatz (|3|4|) 
ll( . which provides the explicit expression of the radial 
distribution function gt{r). By use of this solution, the 
right-hand side of the HRT equation ^ can be written 
in terms of At or, by use of Eq. ([5]), in terms of the free 
energy density At, leading to a closed non linear partial 
differential equation. 

It is apparent the similarity between the Self Consis- 
tent Ornstein Zernike approximation (SCOZA) [T^j and 
this novel formulation of the HRT approach: both theo- 
ries satisfy Eqs. ()3I4I5|) and in both cases use is made of 
the analytic solution of the OZ equation for a Yukawa po- 
tential. Moreover, the consistency condition at the basis 
of SCOZA may be written in the form ([2]) for a specific 
choice of the turning-on procedure of the attractive in- 
teraction: in SCOZA the t parameter (often identified 
with the inverse temperature) only affects the amplitude 
of w(r) without changing its range. This seemingly mi- 
nor difference has profound implications on the behavior 
of the theory in the critical region and close to the phase 
boundary. 

The cut-off function '(/;(<) in Eq. ^ has been chosen 
in order to guarantee the numerical stability of the HRT 
differential equation: il}{t) = (1 + zt/2)^'^ for t < t* 
and -0(0 = (cosh i)~^ for t > t*, where t* is defined by 
imposing the continuity of 4'{t) [l3|. Full details will be 
given in a forthcoming publication. 

Close to the critical point and in the t — > cx) limit, 
the HRT equation ^ simplifies and, when the thermo- 
dynamical variables are properly rescaled, acquires a RG 
structure. The precise form of the rescaled equation ac- 
tually depends on the specific form of the attractive part 
of the interaction (/)(r) ([6|), which in HRT plays the role 
of the smooth cut off of the RG approach. A fixed point 
analysis, similar to that performed for a <I>^ theory 
shows that the HRT equation satisfies scaling and hyper- 
scaling with the non-classical critical exponents shown in 
Table 1. The numerical solution of the full HRT equa- 
tion ^ allows to justify the fixed point analysis on mi- 
croscopic grounds: very close to the critical point, the 
quantity 
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when multiplied by the rcscaling factor e^*, falls into the 
basin of attraction of the fixed point, as shown in Fig. [T] 
However, at very long wavelength (i.e. for t oo), the 
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4.76 
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0.01 


0.332 


1.328 
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4.16 



TABLE I: HRT estimates of the critical exponents and com- 
pressibility amplitude ratio in three dimensions for z = 1.8 
compared to the exact values [T3 ]. 
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FIG. 3: Log-Log plot of the inverse compressibility as a func- 
tion of the reduced temperature above (circles) and below 
the critical temperature along the low-density (triangles) and 
the high-density (squares) branch of the binodal for z = 1.8. 
Solid lines show the expected power law behavior defined by 
the exponent of Table 1 (7 = 1.328). 



FIG. 1; Rescaled inverse compressibility for z = 1.8 during 
the t-evolution at the critical density and reduced tempera- 
tures {T-Tc)/Tc = 1.1 X 10"^2.6 X 10"^-1.4 x 10"^ from 
top to bottom. The dotted line is the fixed point critical value 
and the dashed line is the "zero temperature" fixed point ob- 
tained by a RG analysis of the HRT equation. 



small deviations from the critical point drive the system 
either to the "infinite temperature" or to the "zero tem- 
perature" fixed point. The physical meaning of 
apparent from its definition: it is proportional to the in- 
verse compressibility ([5]) supplemented by the mean field 
contribution associated to the residual part of the poten- 
tial w{r) — Wt{r). In order to understand how the singu- 
larity associated to the first-order liquid- vapor transition 
develops within HRT, it is instructive to follow the "evo- 
lution" of the inverse compressibility Xt~^ foi' ^ fixed tem- 
perature below the critical point and different values of t 
(see Fig. [5]). Long wavelength fluctuations force the in- 
verse compressibility to vanish identically in the t —^ 00 
limit inside the binodal, as customary in the HRT ap- 




FIG. 2: Snapshots of the inverse compressibility along the 
density axis for z = 1.8 and T = 1.1 at three different values 
of the parameter t: t = (dotted line), t = 0.2 (dashed line), 
t = 2.6 (dot-dash line) and t 00 (soUd line). 



proach The novel feature displayed by the smooth 
cut-off formulation, and clearly visible in Fig. [21 is the 
jump of across the phase boundary. A closer inspec- 
tion of the t-evolution of x~^ at coexistence reveals that 
the approach toward zero proceeds differently deep inside 
the binodal, where x~^ remains negative throughout the 
evolution, and close to the phase boundary, where long 
wavelength fluctuations first drive the system towards 
stability (x~^ > 0) and then push the inverse compress- 
ibility to zero. It is tempting to identify the boundary 
between these two regimes as the boundary of metasta- 
bility, i.e. the spinodal curve. A power law fit of the 
compressibility in the critical region x ^ C \T — Td^'^ , 
shown in Fig. [3l is fully consistent with the critical ex- 
ponent 7 reported in Table 1 both above and below the 
critical temperature. The amplitude ratio C+/C_, also 
shown in Table 1, agrees well with the field theoretical 
expectation. The phase diagram of a Yukawa fluid with 
2; = 1.8 and z = 5 is compared with Monte Carlo sim- 
ulations and SCOZA results in Fig. [H The constant 
volume specific heat can be obtained by differentiation 
of the free energy density at convergence Aoo (p, T) . In 
Fig. [5] we plot the specific heat per particle C„ along 
the critical isochore: the data below the critical temper- 
ature are obtained by use of the free energy inside the 
coexistence curve, which is directly accessible in HRT. A 
typical radial distribution function of a z = 1.8 Yukawa 
liquid as obtained by HRT is shown in Fig. [5l The core 
condition is exactly fulfilled due to Eq. Q and the use 
of the analytical solution of the OZ equation. However, 
the linear dependence of c(r) on the attractive interaction 
w(r) implied by the MSA-like closure ^ limits the accu- 
racy of HRT for the description of the local structure of 
the fluid, particularly at low density. Nevertheless, the 
accurate treatment of the physics underlying the first- 
order transition has important consequences on the form 
of g{r) inside the binodal. In the two-phase region, den- 
sity correlations are linear combinations of those of the 
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FIG. 4: Coexistence curve as obtained by HRT (open circles), 
SCOZA (solid line) and Monte Carlo (bars) for z = 5 (upper 
panel) [T^] and z = 1.8 (lower panel) [ig]- The HRT spinodals 
are also shown (full circles). 
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FIG. 5: Left panel: HRT specific heat along the critical 
isochore p = 0.3605 for z = 3. Lower-right panel: radial 
distribution function for z — 1.8, T = 1, p = 0.8. Line: 
HRT result; points: Monte Carlo simulations [l^]. Upper- 
right panel: r oo limit of p^ g{r) inside the coexistence 
curve for z = 1.8 and T = 1. 



two stable phases (liquid L and vapor V). This imphes 
the exact relation linir^oo P^gii^) = p{pv + Pl) — PvPl 
which is satisfied by the numerical solution of the smooth 



cut-off HRT equation, as displayed in Fig. [5l 

In summary, we have shown how a smooth cut off for- 
mulation of HRT provides a consistent picture of the 
equilibrium thermodynamics of a fluid, including the 
complex, singular behavior at first- and second-order 
phase transitions, a distinction between unstable and 
metastable states, and quantitative predictions for the 
coexistence curve, equation of state and specific heat. 
The present formulation of the theory is specific to a sin- 
gle Yukawa interaction, but the available analytical solu- 
tion of the OZ equation for a sum of an arbitrary num- 
ber of Yukawa potentials foreshadows possible generaliza- 
tions. A better representation of correlations may be also 
achieved within the HRT framework, either by adopting 
parametrizations more elaborate than Eqs. p|4p . or by 
closing the hierarchy at the level of the second equation, 
which embodies the effects of density fluctuations on the 
structure of the fluid 6] . Ap plic ations to binary mixtures 
[l^ and non- uniform fluids 17 1 are also possible. 
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